Simulation of interacting fermions with entanglement renormalization 
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We propose an algorithm to simulate interacting fermions on a two dimensional lattice. The ap- 
proach is an extension of the entanglement renormalization technique [Phys. Rev. Lett. 99, 220405 
(2007)] and the related multi-scale entanglement renormalization ansatz. Benchmark calculations 
for free and interacting fermions on lattices ranging from 6 x 6 to 162 x 162 sites with periodic 
boundary conditions confirm the validity of this proposal. 
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The simulation of interacting fermions is of capital 
importance for our understanding of strongly correlated 
phenomena such as high-temperature superconductivity 
or the fractional quantum Hall effect. Quantum Monte 
Carlo techniques, so successful in addressing bosonic 
many-body problems, fail for fermionic systems due to 
the so-called sign problem Many alternative tech- 
niques have been used, including exact diagonalization, 
density matrix renormalization group, dynamical cluster 
approximation, or variational and Gaussian Monte Carlo 
methods But despite all these approaches, even the 
ground state properties of basic lattice models for in- 
teracting fermions, such as the Hubbard model, remain 
highly controversial in two spatial dimensions. 

Recently, entanglement renormalization Q has been 
proposed to efficiently simulate quantum systems on a 
lattice. By means of a coarse-graining transformation, 
the size of the system is progressively reduced until ex- 
act diagonalization becomes feasible. The key of the ap- 
proach is the use of disentanglers to remove short-range 
entanglement at each coarse-graining step. In this way 
the low energy properties of the system are preserved 
while the computational cost is kept under control. An 
approximation to the ground state of the lattice is then 
encoded in the multi-scale entanglement renormalization 
ansatz (MERA) , from which one can compute the ex- 
pected value of local observables and correlators. The 
scheme is scalable and has been used to address arbitrar- 
ily large, two-dimensional spin systems 

In this paper we show that entanglement renormaliza- 
tion can be adapted to address fermionic systems. We 
first explain the key idea underlying the fermionic ver- 
sion of this approach and then demonstrate its validity 
by computing the ground state of fermionic systems on 
lattices of up to 162 x 162 sites. As in the bosonic case, 
the cost of simulations does not depend on whether the 
particles interact, but rather on the amount of entagle- 
ment present in the ground state. 

Fermions on a lattice. — We consider a system of 
fermions on a lattice C made of N sites, where each site 
r 6 C is described by a single fermionic operator c r , with 

{c r -,c\} = S rs I r , {c r ,c s }=0. (1) 

[This simple setting is later generalized.] The system is 



further characterized by a parity preserving Hamiltonian 
H that decomposes as an even polynomial of fermionic 
operators (see e.g. Eqs. lUO, that is as a sum of even 



powers of the c r 's, such as c r c s 



Or ClCrClCs 



The 



occupation number basis for the 2 Ar -dimensional vector 
space of the lattice is given by 

|ii* 2 = (c\YH4T 2 •■• (4)«|00 •■•0), (2) 

where i r = 0,1 indicates the absence/presence of a 
fermion on site r and the vacuum state 00 • • • 0) cor- 
responds to not having any fermion in the lattice. With 
the usual Jordan- Wigner transformation 



= [Z\ ■ ■ ■ Z r -i)o~ r , Z r — I r — 1c\.c 



(3) 



fermionic operators are mapped into spin operators oy, 

[a r ,a' s ] = 6 rs Z r , [a r ,cr s } = 0, (4) 

where Z r acts diagonally on the occupation number ba- 
sis of site r, Z r \0 r ) = |0 r ) and Z r \l r ) — — |l r ), and 
thus preserves the parity of the occupation number of 
a state, whereas the spin operator ay changes parity, 
ay|l r ) = |0 r ), ay|0 r ) = 0. When expressed in spin vari- 
ables, some operators, such as c\c r c\c s = a\a r a\a s , re- 
main supported on just the same sites. Instead, other 
operators, such as cj.c s , develop a string of Z's, 



C„ Ca 



al(Z r+ i ■ ■ ■ Z s -i)a s 



(5) 



In the latter case we say that the bosonic support of the 
operator is larger than its fermionic support. 

Coarse-graining transformation. — Following the 
formalism of entanglement renormalization [H, H|, our 
goal is to coarse-grain the lattice C into a smaller lat- 
tice £ . This is achieved in two steps, as exemplified in 
Fig. Hfi-iv). First, disentanglers u are applied across the 
boundary of blocks of sites of C. Second, isometries w 
are used to map each block of sites of £ into a single site 
of the coarse-grained lattice £'. As extensively discussed 
in Refs. a fundamental aspect of entanglement 

renormalization, on which the efficiency of the approach 
rests, is that local operators remain local under succes- 
sive iterations of the coarse-graining. 
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FIG. 1: (Color online) (i) A 4x4 lattice £ is coarse-grained by 
first (ii) applying disentanglers u on blocks of four sites and 
then (iii) isometries w also on blocks of four sites, producing 
(iv) a 2 x 2 lattice Notice that the Jordan- Wigner order for 
the sites is such that each isometry w acts on four consecutive 
sites. 



Here we shall argue that, with a proper choice of disen- 
tanglers and isometries, locality of operators is preserved 
also in the fermionic case. Two main difficulties need to 
be addressed. On the one hand, fermionic operators such 
as the hoping term c\c s in Eq. [5j have non-local bosonic 
support, with a string of Z's that may involve up to O(N) 
lattice sites. Such non-local supports could in principle 
preclude the use of entanglement renormalization, an ap- 
proach based on transforming local operators only. It 
turns out, however, that this difficulty can be circum- 
vented by focusing on the fermionic support of such op- 
erators, which is local. On the other hand, the locality of 
fermionic supports is not preserved by the type of disen- 
tanglers and isometries used to coarse-grain bosonic sys- 
tems. This can be understood by observing that such dis- 
entanglers and isometries are themselves bosonic (com- 
muting) operators, and bosonic and fermionic (anticom- 
muting) operators are mutually non-local. Fortunatelly, 
this difficulty can be resolved by replacing bosonic disen- 
tanglers u and isometries w with fermionic counterparts, 
that is, with tensors u and w that are local when written 
in fermionic variables. This is further illustrated next 
with an explicit example by considering the simple 4x4 
square lattice of Fig. [TJ 

Fermionic isometries. — Let us first assume that 
the disentanglers in Fig. HJij) arc the identity operator, 
u = I, that is, lattice £ is obtained from C by the use of 
isometries w that map four sites into one. In this case the 
isometries, together with a top tensor, form an ansatz for 
the ground state of H known as tree tensor network 
(TTN) [|], see Fig. Hj). Notice that each isometry w is 
parity preserving (i.e. built as an even polynomial of c 
and c' operators) that act on a block of four consecutive 
sites, e.g. sites 1,2,3,4 e C. Therefore, the fermionic 
and bosonic supports of w coincide (that is, when w is 
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FIG. 2: (Color online) (i) TTN to approximate the ground 
state |\&) of H obtained by adding a top tensor (representing 
the state of £') to the isometries w of Fig. [T] Notice that 
the sites (vertical lines) are arranged in ID according to the 
Jordan- Wigner order. Isometries w are local even when writ- 
ten in terms of spin operators, (ii) MERA obtained by adding 
disentanglers u to the TTN according to Fig. [T] Disentan- 
glers u are delocalized and decompose into sums of terms that 
contain strings of Z's (represented with ribbons). 



expressed in terms of spin operators, there are no strings 
of Z's leaving the block). If we assume (temporarily) 
that the coarse-grained site V £ £' is also described by a 
two-dimensional space with corresponding operator Zy , 
then the parity symmetry of w implies 



w = (ZiZ 2 Z 3 Z 4: )wZ v , 



(6) 



and strings of Z's simply 'commute' with the coarse- 
graining transformation, see Fig. [3]j). 
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FIG. 3: (Color online) (i) A fermionic isometry w preserves 
parity, cf. Eq. [5] (ii) By construction, an isometry w is 
such that the pair w ur anihilates into the identity /. (iii) 
Coarse-graining of the operator c 2 ci2 into Ay Z 2 > By , cf. Eq. 
[71 by using properties (i) and (ii). (iv) The operator C2C12 also 
commutes with a disentangler u if c\cv2. and u have disjoint 
fermionic support. 



It follows that, under coarse-graining, an operator with 
local fermionic support, say c\cyi^ is transformed into an 
operator whose fermionic support is also local, 

a\ (Z 3 Z 4 ■ ■ ■ Zn) (712 -> AyZyBy, (7) 



3 



where Ay and By are parity-changing operators, a prop- 
erty they inherit from <j\ and tJi2, see Fig. [3] (hi). Im- 
portantly, Ay and By are obtained by coarse-graining 
operators a\ and Oyi respectively, whereas the original 
string of Z's simply shrinks into ■ In other words, in 
spite of the presence of a string of Z's, all the manipula- 
tions involved in the coarse-graining of c\c\2 by fermionic 
isometries w can be performed locally and thus efficiently 

Fermionic disentanglers. — Let us now consider 
non-trivial fermionic disentanglers u ^ I, which together 
with the isometries w and a top tensor constitute the 
MERA for the ground state *) of H, see Fig. Wji). 
Again, fermionic disentanglers are built as even polyno- 
mials of c and operators, but since they are not sup- 
ported on consecutive sites of £, they include strings of 
Z's when written in terms of spin variables. 

One may fear that such strings of Z's may turn local 
operators into highly non-local ones. However, this is not 
the case for operators with a local fermionic support. As 
it can be easily checked from Eq. [I] two even polynomi- 
als of c and operators with disjoint fermionic support 
commute with each other. That means, for instance, that 
the operator c^c\2 in Eq. [7] commutes with a disentan- 
gler u with fermionic support on sites 4, 7, 10, 13 6 £, 
m(c2Ci2)m^ = c\c\2, see Fig. That is, fermionic 

disentanglers only expand the fermionic support of oper- 
ators as much as bosonic disentanglers do with bosonic 
supports. In other words, computing an expected value of 
a fermionic operator involves the same number of isome- 
tries/disentanglers as in the bosonic MERA, so that lo- 
cal observables can be computed efficiently. Finally, we 
notice that all the previous considerations still apply in 
lattices where each site is described by a larger vector 
space V, by decomposing the space V r of site r into even 



and odd parity subspaces V r = Yr°^ © YV, with pro- 
jectors Pr°^ and Pr X \ and defining the Z r operator as 

7 = p(0) _ p(l) 

In summary, the use of fermionic disentanglers and 
isometries allows us to maintain the locality of fermionic 
operators during coarse-graining. To put it in the lan- 
guage of quantum circuits, in which the bosonic MERA 
was originally formulated: the causal structure of the 
MERA, consisting of past causal cones with finite 'width' 
|U, is preserved when replacing bosonic wires and gates 
with fermionic ones. As a result, both the TTN [§| and 
MERA [f| @, i| algorithms for two dimensional systems 
can be extended to fermions. Recall that while a TTN 
accurately describes small two dimensional lattices, scal- 
able simulations are only possible with the MERA. See 
Ref. 10] for a thorough technical description of the 
fermionic TTN and MERA algorithms. 

Benchmark calculations. — To test the validity of 
the approach, we first consider a system of free spinless 
fermions with Hamiltonian 
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FIG. 4: (Color online) Left panel: Phase diagram of the free 
fcrmion model JSJ. Right panels: Error in the ground state 
energy obtained with TTN and MERA simulations of a 6 x 6 
lattice with PBC. The lines correspond to different values of 
the refinement parameter \, as indicated in the legend in the 
left panel. The entanglement entropy S1/2 for one half of the 
lattice is larger for those values of A and 7 that lead to larger 
errors. 



on a 6 x 6 lattice with periodic boundary conditions. This 
exactly solvable model exhibits a critical (p-wave) super- 
conducting phase for 7 > 0, < A < 2. and a gapped su- 
perconducting phase for 7 > 0, A > 2 [ll|,LL2]. For 7 = 
the pairing potential vanishes and the model corresponds 
to a free fermion system, i.e. a metal for < A < 2 and 
a band- insulator for A > 2. Fig. [4] shows the error in 
the ground state energy as a function of 7 and A, for 
increasing values of the refinement parameter \, which 
is the dimension of the space V of a coarse-grained site. 
Both TTN and MERA reproduce several significant dig- 
its of the exact solution. |13| The entanglement between 
two halves of the lattice, as measured by the entropy 
S1/2 = — t r (plog 2 p) of the reduced density matrix p for 
half of the lattice, is also plotted. It shows that harder 
computations (those requiring larger values of \ m order 
to achieve a fixed accuracy in the ground state energy) 
correspond to ground states with more entanglement. 
Next we add a nearest neighbor repulsion term to H ttBe , 



#i„t = Hh 



V) clcrclc s 



(9) 



(rs) 



for which an analytical solution no longer exists. We 
emphasize that the algorithm does not require any par- 
ticular modification in order to deal with the interaction. 
Fig. [5]illustrates the convergence of the energy with \ for 
different interaction strengths V, with 7 = 1 and A = 2. 
For small and large interaction we observe a similar con- 
vergence behavior as in the non-interacting case. For an 
interaction strength of the order of the hopping ampli- 
tude, V ~ t = 1, the convergence with x is slower (the 
ground state is again more entangled) but about four 
digits of accuracy seem to still be achieved for large \. 
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FIG. 5: (Color online) Convergence of the ground state energy 
of interacting spinless fermions ([9]) on a 6 x 6 lattice with PBC 
for 7 = 1, A = 2, and varying interaction strength V . The plot 
shows AE = E x — -E X max) the difference between the energy 
as a function of \ (squares for TTN, circles for MERA) and 
our best MERA result E Xsalac , where Xmax = 120. Again, the 
entanglement entropy S1/2 shows a strong correlation between 
ground state entanglement and convergence in \. 
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FIG. 6: (Color online) The error in the ground state energy as 
a function of the linear size L of a square lattice, obtained with 
the fermionic MERA algorithm with x = 4 (non-interacting 
case, 7 = 1), is of the same order of magnitude for small and 
large systems, even in the critical regime A < 2. Simulations 
with interacting fermions (not plotted) show an analogous 
pattern of energies, but there are no exact results to compare 
with. 



Finally, Fig. [5] shows the error in the ground state en- 
ergy for lattices of up to 162 x 162 sites, demonstrating 
the scalability of the fermionic MERA algorithm in two 
spatial dimensions. 



Discussion. — We have shown that interacting 
fermionic systems can be addressed within the formal- 
ism of entanglement renormalization. Importantly, as 
with spin systems, the cost of simulations is not deter- 
mined by the strength of interactions, but by the amount 
of entanglement in the ground state. While a precise 
characterization of two-dimensional fermionic systems in 
terms of ground state entanglement is still missing, our 
results for interacting fermions are consistent with those 
obtained in Ref. for free fermions and suggest that, 
broadly speaking, gapped systems are the easiest to sim- 
ulate. These are followed by gapless systems with a finite 
number of gapless modes, such as the critical supercon- 
ducting phase in Fig. Gapless systems with a ID 
Fermi surface, the most entangled systems, appear also 
as the most challenging. 

We envisage that the present approach will help ad- 
dress long standing questions in strongly correlated sys- 
tems. Presently, a major limitation is due to the scaling 
0(x 16 log L) of the computational cost (translation in- 
variant case [EEl)- While the mild dependence on L 
ensures scalability, only small values of x can be consid- 
ered. For instance, each L = 162 simulation of Fig. [B] 
with x — 4 already took several days on a 3GHz dual- 
core desktop PC with 2Gb of RAM [HI. A number of 
strategies are being considered to improve the computa- 
tional cost, such as alternative coarse-graining schemes, 
exploitation of internal symmetries (e.g. particle conser- 
vation or spin isotropy), use of parallelized code on larger 
computers, and variational Monte Carlo sampling tech- 
niques. Work in progress includes exploring the ground 
state phase diagram of the Hubbard model. 
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cussions, and S. Haas, L. Ding and N. Ali for clarifi- 
cations concerning the free fermion model ([5]). Support 
from the Australian Research Council (APA, FF0668731, 
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Note added. Short after this work was made avail- 
able online, a largely equivalent approach was indepen- 
dently presented by C. Pineda, T. Barthel and J. Eisert 
in Ref. 0. 
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